Extended Hartree-Fock method based on pair density functional theory 
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A practical electronic structure method in which a two-body functional is the fundamental variable 
is constructed. The basic formalism of our method is equivalent to Hartree-Fock density matrix 
functional theory [M. Levy in Density Matrices and Density Functionals, Ed. R. Erdahl and V. H. 
Smith Jr., D. Reidel, (1987)]. The implementation of the method consists of solving Hartree- 
Fock equations and using the resulting orbitals to calculate two-body corrections to account for 
correlation. The correction terms are constructed so that the energy of the system in the absence 
of external potentials can be made to correspond to approximate expressions for the energy of 
the homogeneous electron gas. In this work the approximate expressions we use are based on 
the high-density limit of the homogeneous electron gas. Self-interaction is excluded from the two- 
body functional itself. It is shown that our pair density based functional does not suffer from the 
divergence present in many density functionals when homogeneous scaling is applied. Calculations 
based on our pair density functional lead to quantitative results for the correlation energies of atomic 
test cases. 
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I. INTRODUCTION 



The most commonly used methods in electronic struc- 
ture for atoms and molecules are density-functional the- 
ory P, U, i, H (DFT) and Hartree-Fock [s', ^ (HP) the- 
ory. The former is based on the Hohenberg-Kohn the- 
orems, which state that all ground state observables, 
in particular the energy can be written as a functional 
of the one-body density [l|. Since the functional itself 
is unknown, in actual implementations the system of 
interacting electrons is mapped onto an auxiliary sys- 
tem of non-interacting electrons resulting in the Kohn- 
Sham single-particle equations [2]. This ansatz implic- 
itly assumes that all densities originating from anti- 
symmetric many-body wave-functions can be represented 
by a non-interacting wave- function. The Kohn-Sham 
energy functional includes a correction term, known as 
the exchange-correlation energy, which accounts for ef- 
fects of exchange and correlation, and is usually approx- 
imated usingthe theory of the homogeneous electron 
gas i, 0, BM Ei [O, [ll [11 (HEG). This correction 
term usually depends on the one-body density, although 
orbital-dependent variations also exist |14| . 

The HF approximation [5, 6] is based on the use of 
a Slater determinant for the many-body wave-function, 
hence the work-horse equations are also single-particle 
equations (as in DFT), but exchange is incorporated in 
this case. Correlation effects are neglected. There are a 
variety of methods which incorporate correlation effects 
into HF, such as use of a linear combination of Slater de- 
terminants [6| , perturbation theory [a, '6| , or a combina- 
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tion thereof, but all at the cost of greatly increasing com- 
putational demand 15] (both methods involve the use of 
virtual orbitals, increasing the size of the Hilbert space in 
which the wave-function is calculated). The connections 
between HF and DFT have also been investigated [16!, llll 
and it has been shown that the Hohenberg-Kohn theo- 
rems apply within the Hilbert space of Slater determi- 
nants, hence it is in principle possible to construct a DFT 
that corresponds to HF. 

There has also been active interest in extending DFT 
by usingthepair densityjlk 19, 20, 21, 22, 23, 23, [H, 
Hi, [13, M [Si, [13, [M Hill] (pair density functional 
theory (PDFT)) or the density matrix [H, [ll, [H, [13, [11] 
as the fundamental variable. Most PDFT studies have 
dealt with formal issues, an ex amp le of an exception be- 
ing the work of Gonis et al. |30| where an implemen- 
tation is developed in the context of strongly correlated 
systems. The development of a practical electronic struc- 
ture method based on PDFT, and applicable to atomic 
and molecular systems in general appears to be lacking. 

PDFT is an extension of the standard DFT in that 
the Hohenberg-Kohn theorems \M a re g eneralized to the 
two-body (or TV-body) density [ll, \m, [13, HH [llpl, 
[2^, [23, [so] ■ Application of the Kohn-Sham ansatz [^ is 
not straightforward in the case of PDFT, since a non- 
interacting pair-density, such as that obtained from a 
Slater determinant can not represent all possible pair 
densities arising from an arbitrary anti-symmetric wave- 
function. The lack of correlation between electrons with 
anti-parallel spin presents an obstacle to representing an 
arbitrary anti-symmetric many-body wave-function via 
Slater determinants. 

An approach related to PDFT based on reduced den- 
sity matrices (RDM) is also under considerable investiga- 
tion [Hm, [11,133, [si- Levy has demonstrated [H that 
a correction term due to correlation is a unique functional 



of the Hartree-Fock density matrix. Levy has also sug- 
gested [331 possible starting points for the construction 
of such correlation corrections. 



II. FORMAL DEVELOPMENT 



A. Two-body density-functionals 



In this work we construct a PDFT-based method that 
is robust, simple, and it can be implemented in an HF 
context. We argue that an exact energy functional can 
in principle be constructed based on the HF pair den- 
sity. Since the HF pair density and one-particle RDM 
are related one-to-one, our formalism is equivalent to 
that derived by Levy. [35| A pair density based correction 
to HF has also recently been suggested by Higuchi and 
Higuchi [33. 



We also construct approximations to the exact cor- 
relation functional via generalizing known scaling rela- 
tions [l9|, which allow the construction of pair density 
dependent correlation functionals based on the theory 
of the homogeneous electron gas (HEG) 0, 0, B S [l3, 
[ill . [l3 . [l3l . The method suggested by Higuchi and 
Higuchi j3] uses the standard scaling relations derived 
by Levy and Ziesche [19|. In particular, in this work, 
we construct correlation functionals which correspond 
to expressions valid for the HEG in the high-density 
limit 1^, [13, [nl, [l3, [3 (expected to be a good approxi- 
mation for the atomic test cases presented herein). Our 
procedure can, in principle, be extended to other func- 
tional forms (power series in the density, such as the low- 
density form) for the correlation energy es well. Since our 
method amounts to implementing a correction term to an 
HF calculation (to account for correlation), and it does 
not require the use of virtual orbitals, the scaling with 
system size is that of HF itself. We apply our method to 
calculate atomic energies and obtain quantitative agree- 
ment with known correlation energies. 



We also investigate the properties of our constructed 
functionals under homogeneous coordinate scaling. DFT 
correlation energy functionals with logarithmic terms [9, 
llO . [ill . [l2 . [l3l | are known to be divergent when the scal- 
ing parameter is made infinit e [ 39l . [40|. This behavior 
contradicts the known [40, [4l|, \^ scaling behavior of the 
actual correlation energy which is bounded below. It is 
shown that our correlation functionals do not exhibit this 
deficiency. 



This paper is organized as follows. In Section[TT|the for- 
mal development is presented. We describe a procedure 
to construct correlation functionals based on known scal- 
ing relations in PDFT, and comment on the properties of 
the kinetic energy functional. The explicit construction 
of energy functionals is explained in section IIIII In sec- 
tion llVl we discuss the details of our calculations. Section 
[V| consists of our results and analysis. We conclude our 
work in section FVll 



The ground state energy of a many-body fermionic sys- 
tem has been shown to be a unique functional of the 
single-particle many-body density matrix [j, [3J] . In Ref. 
I35I it was shown that there exists an additive correction 
to the HF energy functional that depends on the HF 
one-body RDM. The proof is based on investigating the 
correspondence between the HF one-body RDM and the 
one-body potential of a given system. In particular an 
exact formula is provided by which one can determine 
the external potential as a function of the HF orbitals. 

In standard DFT, an applicable method is constructed 
via the Kohn-Sham ansatz [2|. This ansatz involves 
a mapping between the interacting system of electrons 
to an auxiliary non-interacting one, enabling the use of 
single-particle equations. The solutions of the single- 
particle equations are single-particle orbitals, and one 
must choose a way of writing the density in terms of 
these single particle orbitals. The most often used choice 
is the one that satisfies the Pauli principle 



p(i') = Xl-^'^'^i'?^'^'^w 



(1) 



where p(r) denotes the density fi^a- denotes an occupa- 
tion number and (j)i^cr (r) denotes a particular spin-orbital. 
Usually fi^a = 1 which corresponds to the unrestricted 
HF density. To our knowledge it is not proved that the 
non-interacting form for the density based on a single 
determinant can represent all interacting densities corre- 
sponding to anti-symmetric many-body wave- functions. 
While standard DFT was made into an applicable elec- 
tronic structure method by the Kohn-Sham ansatz, the 
generalization of the ansatz to the case of two-body func- 
tionals such as the pair density is not strai ghtf orward, 
due to representability issues [H [13, HI, (H,!!!. Writ- 
ing a non-interacting energy- functional requires the defi- 
nition of the pair density n(r, r') in terms of the orbitals 
of the non-interacting system. An obvious initial choice 
is 

nHF{r,r') = ^ {|0,.(r)n0,,<,,(r')P 

- 0,,.(r)0,,,(r')(/.j,.'(r)0,-,,(r')} (2) 

obtained from a Slater determinant (also known as the 
Hartree-Fock pair density). When this definition for the 
pair density is used, it can be shown that the motion 
of electrons with anti-parallel spins is uncorrelated ^5|]. 
Since this is not so for the anti-symmetric many-body 
wave- function, exact representation of an interacting sys- 
tem by a non-interacting one does not appear to be possi- 
ble. In the theory of correlated wave- functions this effect 
is usually implemented via the cusp condition for corre- 
lation factors [§|. 





Gell-Mann and Brueckner 


Nozieres and Pines 


RPA 


A, 


-0.048 


-0.058 


-0.071 


Be 


+0.0311 


+0.016 


+0.0311 


A 


-0.0619 


-0.0649 


-0.0859 


B 


-0.0104 


-0.00517 


-0.0104 


A 


-0.0598 


-0.0639 


-0.0838 


B 


-0.0622 


-0.031 


-0.0622 



TABLE I: Table of constants for different correlation energy 
functionals (Har) (functional forms are given in Eqs. (|20|l and 



On the other hand, using the result of Levy [35| one 
can easily argue that the energy is a unique functional 
of the Hartree-Fock pair density since the density matrix 
and the Hartree-Fock pair density are in a one-to-one 
relation. Hence the correlation energy may be written 



Ec[nHF] = E[nHF] - EnFlnHp], 



(3) 



i.e. as a unique functional of the HF two-body den- 
sity (in Eq. ([3]) E[nHF][EHF[nHF]) denote the ex- 
act energy(Hartree-Fock energy) as a functional of the 
Hartree-Fock pair density n\^HF])- Below we construct 
approximations for Ec[nF[F]- 

Here, as in the case of DFT, arguments based on map- 
pings establish one to one relations between observables 
(ground state energy) and coordinate-dependent func- 
tionals (density), however the exact functional relations 
are not given by the formalism. While in the Kohn- 
Sham DFT method a mapping between a non-interacting 
and interacting system is invoked, where the two are 
constrained to have the same density, in the formalism 
developed above one does not have the two-body den- 
sity readily available. Although, in principle, a pertur- 
bative approach carried out to infinite order gives an 
exact solution, in practice one has to content oneself 
with an approximate solution to this problem. Higuchi 
and Higuchi [33| suggest minimizing the energy of an 
augmented Hartree-Fock equation, but this method can 
only be approximate due to lack of opposite-spin corre- 
lation. Possible alternative approaches are perturbation 
theory [ifl, llllj E^] j two-body generalization of plane- wave 
perturbation theory [ll|, or variational approaches such 
as the Fermi hyper- netted-chain [4J, |45| method. A two- 
body density can also be obtained from an HF solution 
via constructing a Jastrow-Slater wave-function to ac- 
count generally for correlation, and the relevant observ- 
ables could be evaluated by quantum Monte Carlo meth- 
ods [4y|. In all of these methods, the HF wave- function 
serves as input, and the output consists of the many-body 
properties. 

Unfortunately our formalism does not provide an in- 
version formula that can be used to determine the two- 
body potential from the two-body density in the HF case. 
This is due to the fact that, although the HF two-body 



density is a two-body object, it is composed of single- 
particle orbitals. If a particular two-body potential is as- 
sumed (Coulomb repulsion between electrons) then the 



35| can be used to 



inversion formula presented by Levy 
obtain the one-body potential. 

The method constructed from the above reasoning re- 
quires an HF calculation, supplemented by a correction 
term due to correlation. The correction term is a func- 
tional of the HF pair density. Knowledge of the HF wave- 
function gives exact information about the ground state 
energy (and by extension all ground state observables) in 
principle, and it is not necessary to calculate the exact 
wave-function, if one is interested only in ground state 
expectation values. This can not be done exactly, due 
to the lack of the exact functional form, and approxima- 
tion schemes need to be conceived for the particular set 
of observables in mind. In this study we develop such a 
scheme for the ground state energy, or more specifically 
the correlation energy. 



B. Generalized scaling relations and sum rules 

To construct correlation functionals -E'ci^gf] we write 
the scaling relations due to Levy and Ziesche [l^ . In Ref. 
[l9J it was shown that possible terms that contribute to an 
expansion of the kinetic energy functional (per particle) 
of the form 



t\n\ 



N ^ 

b 



Ah / drdr'n(r, r')°|r 



„'ifc 



(4) 



must be such that the relation 6a — 6 = 8 is satisfied. To 
derive the HF kinetic energy functional for the homoge- 
neous case we substitute the relation 



nHF{r,r') = p^g{\r - r'\), 



(5) 



into Eq. ^. In Eq. ([5]) p denotes the one-body density, 
and g{r) denotes the radial distribution function. Us- 
ing the dimensionless distance scaled by the Fermi wave- 
vector X = fci?r, where fci? = (37r^/5)3 , elementary manip- 
ulations lead to the kinetic energy expression 



t[nHF] = p'- 



b 



dxx'+'gixY, 



where 



g{x) = {1 



2x6 



(sin(a;) — xcos(a;))^}. 



(6) 



(7) 



The " on A^ indicates that some constants have been ab- 
sorbed into Ab. It is well-known 0,|j, 0] that the kinetic 
energy in the HF approximation can be written as 



t[nHF] 



10 



(3^2)3p3^ 



(8) 



thus, from Eqs. ([6]) and ([8]) a sum rule is established on 
the coefficients Ai, for the homogeneous non-interacting 
case, which we write as 



E^' 



dxx''+'gixr^^{3n')l 



(9) 



Regarding an homogeneous interacting system we can 
also conclude from the above analysis that if the kinetic 
energy exhibits scaling other than pa then the scaled ra- 
dial distribution function g{x) also depends on the den- 
sity. 

It is possible and useful to generalize the scaling rela- 
tions of Levy and Ziesche [l9| so that a particular func- 
tional can reproduce a given density dependence when a 
homogeneous non-interacting system is invoked. A gen- 
eral functional of the pair density n may be written 



/(r)N = ^ E ^« / '^'■'^'^'^(^' ^')"i^ 



„'|b 



(10) 



where the condition 6a — 6 = 6 -|- F has to be satisfied. 
For the homogeneous non-interacting case (i.e. when Eq. 
([5]) is substituted for uhf) /(r)[?T.] will exhibit a density 
scaling 



f{T)[nHF]^Crp^^ 
resulting in the generalized sum rule 



PCX) 



9[x) 



(11) 



(12) 



valid for the auxiliary non-interacting system. The ~ on 
C{a) indicates that some constants have been absorbed 
into C{a). Using Eq. (fTTj) one can construct pair den- 
sity analogs for correlation functionals that are given in 
terms of power series in the density p (examples are the 
low and high-density approximations). It is also possible 
to obtain sum rules of the form Eq. (fT2|l valid for the 
correlation energy. 



C. Properties of the kinetic energy functional 

The properties of the integral in Eq. ^ can now be 
studied since the functional form of the radial distribu- 
tion function is known (Eq. ^). The behavior of g{x) 
at zero and infinity place bounds on the exponent a (and 
hence b). At large distances 5(x)° approaches 1, hence 



6 + 2< -1. 



(13) 



Close to zero, however g{x) for the unpolarized case ap- 
proaches a constant (i) which leads to the relation 



b + 2>-l, 



(14) 



which is inconsistent with Eq. (J13p . Hence for the spin- 
unpolarized case the integral is divergent. This diver- 
gence is due to the radial distribution function being fi- 
nite at the origin, in the case of the ideal Fermi gas. 



The divergence can be eliminated in several ways. One 
is to assume that the kinetic energy is a sum of terms, 
one for each electron spin component. The radial dis- 
tribution function corresponding to one spin component 
only is zero at the origin, therefore the divergence docs 
not arise. In HF this procedure is implicit. Corrections 
to the energy functional can also be constructed in terms 
of the exchange hole, as is done below. 



III. CONSTRUCTION OF ENERGY 
FUNCTIONALS 

Since the HF method is now well-established, our prin- 
cipal aim is to construct energy functionals that are 
tractable by equivalent numerical methods. To achieve 
this aim, we construct approximations for the correla- 
tion energy Ec in terms of the Hartree-Fock pair density. 
Since in the HF theory of the HEG the quantity that ap- 
pears is the so called exchange-correlation hole defined 
as 



= (r,r')=n(r,r')-p(r)p(r'). 



(15) 



we find it more convenient to use n^c as our fundamental 
quantity (input function). Below we derive approxima- 
tions from the theory of the HEG, and in HEG the energy 
would diverge without the constant positive background. 
The second term in Eq. pS]) corresponds to the positive 
background, hence the divergence is canceled. The scal- 
ing relations derived in the previous section arc valid for 
rixc as well, since both terms on the right-hand side of Eq. 
(fT5)) exhibit the same coordinate scaling. When the HF 
approximation is invoked the exchange correlation hole 
becomes the exchange hole 



nr,{r,r') = 



(l)i,a (r)0i,^ (i"')'/'i:<T (i")'/'i.<T (!■')• (16) 



Eq. P^ reintroduces self-interaction which can be ex- 
cluded here by restricting the summation to i ^ j, 



n^(r,r') = - ^ '/',:,<T(r) 



.(r')^,..(r)^,-,(r'), (17) 



i^i," 



where the prime indicates the self-interaction corrected 
(SIC) sum. In the two-body density functional the- 
ory presented here self-interaction is thus excluded from 
the density functional itself. In DFT SIC is usually 
achieved by subtracting energies (LSD energy function- 
als for Hartree, exchange and correlation) corresponding 
to single-particle densities [ij]. Recently, Lundin and 
Eriksson [43] have proposed a SIC procedure where the 
self-interaction is addressed by subtracting single-particle 
densities from the total density. 

It has been pointed out based on scaling arguments 
in the case of the exchange energy that self-interaction 
terms disappear in the thermodynamic limit [48l . |49j. 
The reason for the disappearance of the self-interaction 



Ec 


He 


Be 


Ne 


Mg 


Ar 


RPA 


0.130 


0.251 


0.884 


1.17 


1.69 


GB 


0.0821 


0.154 


0.644 


0.782 


1.26 


NP 


0.109 


0.213 


0.662 


0.935 


1.24 


RPA - SIC ~ LDA 


0.014 


0.033 


0.198 


0.222 


0.398 


GB - SIC - LDA 


0.014 


0.033 


0.198 


0.222 


0.398 


NP - SIC - LDA 


0.0072 


0.016 


0.0987 


0.111 


0.198 


RPA - SIC ~ LSD 


0.082 


0.161 


0.589 


0.689 


1.13 


GB - SIC - LSD 


0.058 


0.113 


0.469 


0.545 


0.912 


NP - SIC - LSD 


0.063 


0.124 


0.404 


0.476 


0.761 


VMC [57] 


- 


0.073 


0.346 


0.372 


0.576 


Clementi and 
Hoffmann [56.1 


0.045 


0.094 


0.3870 


0.438 


0.722 



TABLE II: Correlation energies (Har) from DFT calculations 
using various correlation functionals and their self-interaction 
corrected versions for closed-shell atoms. 



Ec 


Ca 


Zn 


Kr 


Xe 


RPA 


1.85 


3.15 


3.88 


6.12 


GB 


1.37 


2.43 


3.01 


4.83 


NP 


1.36 


2.23 


2.73 


4.24 


RPA - SIC - LDA 


0.429 


0.822 


1.01 


1.66 


GB - SIC - LDA 


0.429 


0.822 


1.01 


1.66 


NP - SIC - LDA 


0.214 


0.410 


0.505 


0.83 


RPA ~ SIC ~ LSD 


1.23 


2.13 


2.62 


4.15 


GB - SIC - LSD 


0.994 


1.77 


2.18 


3.50 


NP - SIC - LSD 


0.836 


1.39 


1.70 


2.67 


VMC [57] 


0.619 


1.27 


1.43 


2.03 


Clementi and 
Hoffmann [56] 


0.842 


1.74 


2.26 


4.04 



TABLE HI: Correlation energies (Har) from DFT calculations 
using various correlation functionals and their self-interaction 
corrected versions for closed-shell atoms. 



terms is that self-interaction terms have to scale as A'' 
with the number of particles, whereas the overall num- 
ber of interactions between different particles scales as 
N{N — l)/2. An underlying assumption of the argument 
which is vaUd for the exchange energy is a single-particle 
picture in which one can decompose the interactions be- 
tween different pairs of particles. Since the correlation 
energy in the high-density limit which serves as the basis 
for our approximations here is also a single-particle based 
expression (the diagrammatic summation gives an oper- 
ator whose expectation value has to be evaluated over 
a Slater determinant [sl, [l2, Il3|), and it is also valid in 
the thermodynamic limit the arguments in Refs. [43, H 
are valid for the correlation functionals used here (below 
Eq. (fT8|) ) as well. Thus the expressions used below for 
the correlation energy (Eq. (jlSp). valid in the thermody- 
namic limit, do not exhibit self-interaction. It can also be 
expected that use of the SIC exchange hole to reproduce 
the functionals such as Eq. (fTS]) is thus advantageous, 
since it lacks self-interaction by construction. 



edp) = A + Blnp. 



(20) 



In addition to power series forms, it is also possible 
to represent logarithmic terms in a pair density context. 
A simple functional form for the correlation energy per 
particle common to the approximations of Gell-Mann and 
Brueckner [Hi (GB), Nozieres and Pines ^ (NP), and 
the random-phase approximation 9] (RPA) may be writ- 
ten as 



ec{rs) = Ac + BJnrg, 
or taking advantage of the definition of r^. 



47rr^ 



1 
P 



(18) 



(19) 



The first term (constant A) converts analogously to Eqs. 
(HI) and ^. In the local-density approximation for DFT 
without SIC the constant term would correspond to AN, 
where N denotes the number of particles. The second 
term, however, can be written in a PDFT form as 



Blnp = G + ^ / rfrdr'n,(r, r')hi|r - r'| 



(21) 



and hence one can arrive at a PDFT analog of the form 
edfi^] ='^+— drdr'n^(r,r')ln|r-r'|. (22) 

The derivation of the constants proceeds analogously 
to the derivation in Eqs. ([4]) and ([6]) (one uses the HF 
pair density, and the definition of the Fermi wave- vector) . 
The values of the constants for GB and NP for both cases 
(DFT and PDFT) are given in Table H The values of 
the constant A and B are fixed by requiring that the 
correlation energies be equal to the known values for the 
homogeneous case (i.e. Eq. (122]) reproduces Eq. ([20)) ). 
The derivation of the relation between the constants is 
given in the Appendix. 

The use of the SIC exchange-correlation hole in Eq. 
([22[) has two important consequences. It has been 
pointed out [3^, ^^ that under homogeneous coordinate 
scaling the logarithmic term in the approximate correla- 
tion energy functionals such as GB, NP, and RPA, di- 
verges as the scaling parameter A is made infinite. This 
divergence is avoided in Eq. ([22|l due to the use of 
the exchange hole. Indeed, coordinate scaling shifts the 
functional in Eq. (j22p by an amount proportional to 
/ (ir(ir'ria;(r, r')lnA, which integrates to zero in the case 
of an orthonormal basis set. This is advantageous since 
it has also been shown that the actual correlation func- 
tional scales to a bounded constant Ho, Bj], B3l. For 



two-electron atoms the constant has been calculated: - 
0.0467Har [50, [Sll- In our case this constant is zero for 
one class of our constructed functional (Eq. ([23| ) and 
it corresponds to the constant A shown in Table U for the 
other (Eq. (|24p). The divergence would arise had we not 
applied the SIC to the exchange hole. 

Another consequence is that self-interaction is avoided 
by construction. The first term in Eq. (P^ is therefore 
zero, and the constant term A does not contribute. This 
can be shown by considering Eqs. (fTO|) and (fTTj) . When 
the exponent F in Eq. pT|) is taken to be zero, and the 
exponent a in Eq. pop is taken to be one it follows that 
6 = 0. In this case a constant term in the correlation 
functional again corresponds to an integration over the 
exchange hole (Eq. (dH])) giving a contribution of zero. 
The disappearance of the constant term is also present 
in the SIC of Lundin and Eriksson [43 • In the stan- 
dard DFT version of the correlation functional the con- 
stant term also does not contribute if the SIC procedure 
of Perdew and Zunger [ij] is applied to the functional 
within the local density approximation (for an example 
of an application of the SIC in an LDA context see Ref. 
|52| ). since the constant term in the correlation energy of 
the electron gas (which is defined as energy per parti- 
cle) corresponds to a term consisting of an integral over 
the density to the first power. Applying the SIC to such 
an integral leads to cancellation. The problem does not 
arise when the SIC is formulated in the local-spin density 
(LSD) approximation, since the constant is scaled [53] in 
the SIC term which is subtracted from the uncorrected 
functional. Below we also compare using the two differ- 
ent versions of SIC for the constant term. 

It should be emphasized that in both cases the can- 
cellation of the constant term is an artifact of the self- 
interaction procedure used, which, in this case estimates 
the self-interaction component to be as large as the ex- 
change constant itself. The constant in Eq. (fTSl) . A^ 
arises from direct and exchange interactions [12| . and 
it would be a severe approximation to discard it. One 
possible way around this problem in our formalism, in 
principle, is to relax the condition a — 1. Such a pro- 
cedure would, however, significantly worsen the scaling 
of the method, since explicit evaluation of the two-body 
density would be necessitated. Instead we apply the stan- 
dard LSD based SIC procedure of Perdew and Zunger [iJI 
to A in the correlation functionals, which amounts to an 
additive term of AN/2. Applying an LSD based SIC pro- 
cedure to the constant term is not inconsistent with the 
SIC procedure we use, since our SIC corrected exchange 
hole hx is a sum over different spin-components. 

Another point to emphasize is that although the ex- 
change hole includes correlations between electrons with 
parallel spins only, the energy functionals constructed ac- 
cording to our procedure account, at least in principle, 
for all interactions that are included in the approximate 
high-density correlation energy functionals (RPA, GB, 
NP). The interactions accounted for are determined by 
the functionals to which our PDFT based functionals are 



made to correspond. By analogous logic, in DFT the cor- 
relation energy functional depends on the density only, a 
one-body object without correlations, yet it accounts for 
correlation by virtue of being required to reproduce the 
energetics of approximations to correlation in the HEG. 
Our correlation energy functional (total as opposed to 
per particle) thus looks like 



E. 



^'h. 



+ - 



B 



drdr'nx{r, r')ln|r — r' 



(23) 



When the constant term A is subject to LSD based SIC 
treatment the correlation energy functional becomes 



E'Jin.] = ^ + f / drdr'hAr, r')ln|r - r'| 



(24) 



IV. CALCULATION DETAILS 

While the above formalism amounts to a modified HE, 
here we perform a normal HE calculation and calculate 
the correlation energies using the HF orbitals. Thus our 
calculations amount to a perturbative treatment. This 
treatment was also suggested by Levy [Sa]- Given that 
for our test cases (closed sub-shell atoms) the correlation 
energies are ca. 1% of the total energy, the perturbative 
approach can be expected to be a good approximation. 

We have performed atomic HF calculations using the 
most recent version |54| of the atomic program written 
by Fischer |55j . Subsequently we have used the orbitals 
obtained from the atomic program to calculate the corre- 
lation energy contributions given by the PDFT analogs 
that we constructed. For comparison we have also cal- 
culated the correlation energies from the DFT analogs of 
GB and NP. The GB [Ti| and NP [lS\ correlation func- 
tionals are based on a diagrammatic perturbation sum- 
mation for the HEG. While the GB is an expansion valid 
in the limit of high-density, the NP functional is an in- 
terpolation formula approximately valid for all density 
ranges [1, [1]. 

For reference purposes we have calculated the correla- 
tion energies within DFT. In TableslTTIand lllll we present 
the correlation energies for a sequence of closed-shell and 
closed sub-shell atoms calculated using the GB, NP, and 
RPA forms. The results of Clementi and Hoffmann [5y| 
and the variational Monte Carlo results of Buendia et 
al. J57| are also shown. As expected the correlation ener- 
gies are overestimated due to the absence of SIC. Appli- 
cation of SIC based on LDA leads to an underestimation 
of the correlation energy, whereas the LSD based SIC 
leads to an improvement over both previous methods. 
When improved correlation energy functionals are used 
the correlation energy is quantitatively recovered [lj| . 



RESULTS AND ANALYSIS 



E, 


He 


Be 


Ne 


Mg 


Ar 


RPA-I 


0.0 


0.00133 


0.0704 


0.0755 


0.156 


GB-I 


0.0 


0.00133 


0.0704 


0.0755 


0.156 


NP-I 


0.0 


0.000661 


0.0351 


0.0376 


0.0776 


RPA - II 


0.0838 


0.169 


0.489 


0.578 


0.910 


GB~II 


0.0598 


0.121 


0.369 


0.434 


0.694 


NP-II 


0.0639 


0.128 


0.354 


0.421 


0.652 


VMC [57] 


- 


0.073 


0.346 


0.372 


0.576 


Clementi and 
Hoffmann [56,1 


0.045 


0.094 


0.3870 


0.438 


0.722 



TABLE IV: Correlation energies (Har) from PDFT calcula- 
tions using various exchange correlation functionals. For the 
meaning of the different functionals (I, H) see the text. 



E, 


Ca 


Zn 


Kr 


Xe 


RPA-I 


0.166 


0.362 


0.340 


0.829 


GB-I 


0.166 


0.362 


0.340 


0.829 


NP~I 


0.0830 


0.181 


0.170 


0.368 


RPA - II 


0.99 


1.62 


1.84 


3.09 


GB-II 


0.763 


1.26 


1.42 


2.44 


NP-II 


0.722 


1.14 


1.32 


2.13 


VMC [52 


0.619 


1.27 


1.43 


2.03 


Clementi and 
Hoffmann [56] 


0.842 


1.74 


2.26 


4.04 



TABLE V: Correlation energies (Har) from PDFT calcula- 
tions using various exchange correlation functionals. For the 
meaning of the different functionals (I, II) see the text. 



In Table IIVI the correlation energies based on the 
PDFT analogs of the RPA, GB, and NP correlation func- 
tionals are presented for a sequence of closed sub-shell 
atoms He, Be, Ne, Mg, Ar (I refers to the functional 
given in Eq. P5)) ). The calculated correlation energies 
capture the qualitative trend for all three functionals. An 
apparent deficiency is the fact that the correlation en- 
ergy for the helium atom is zero. This deficiency stems 
from the use of the SIC exchange-correlation hole as the 
functional to construct the correlation energy functional. 
The SIC exchange correlation hole (Eq. ^T5\i ) for a non- 
interacting system consists of interactions between elec- 
trons with parallel spins. Since the two electrons in a 
helium atom (ground state) are an anti-parallel pair it 
can be expected that a correlation energy functional of 
the exchange correlation hole give zero as a result. 

We note that there is a qualitative similarity and at the 
same time a quantitative difference between this short- 
coming of functional I (Eq. (^5]) ) and the spurious self- 
interaction in the case of the original LDA density func- 
tionals. In LDA it is a property of the input function 
(the one-body density) that gives rise to the spurious 
contribution. Since the one-body density is some finite 



positive-definite function even for a single electron a func- 
tional that depends on it gives a finite contribution in 
general. The similar artifact of the SIC exchange hole as 
an input function is that it lacks terms with anti-parallel 
interactions, therefore it is not expected to give contribu- 
tions when a system with one single interaction between 
electrons of anti-parallel spins is considered, such as he- 
lium. 

In Table lTVl results for the approximately corrected ver- 
sion (Eq. ([M]) ) of the PDFT correlation energy are also 
presented (RPA-II, GB-II, and NP-II). All three func- 
tionals are improved considerably compared to their val- 
ues without the SIC Table [TVl and quantitative (near ex- 
act) agreement is observed for the GB-II and NP-II func- 
tionals when compared to the results of Clementi and 
Hoffmann ^56] . 

The pair density ansatz GB-II appears to give the best 
results for the sequence of atoms when compared to the 
results of Clementi and Hoffmann [5y|. The results for 
GB-II are also consistently better than the corresponding 
DFT results in Table |TI](GB-SIC-LSD). 

The recent results of Buendia et al. [53| are also pre- 
sented for comparison. This study is based on a varia- 
tional Monte Carlo calculation using a correlated basis 
function. Our results compare well with the known re- 
sults, and overall the agreement appears to be as good 
as the VMC calculation of Buendia et al. [53| . The fact 
that in our case the agreement worsens for small atom 
systems can be attributed to the fact that the correla- 
tion functional is approximated from the HEG, a system 
in the thermodynamic limit. 

In Table |V] results for a heavier sequence of closed sub- 
shell atoms Ca, Zn, Kr, Xe are also presented. The cor- 
relation energies of Clementi and Hoffmann [5y| are also 
shown for comparison. While our results are no longer 
in as good agreement for the correlation energy as for 
the lighter atoms in Table [TVl they appear to be close to 
the results of the VMC calculation of Buendia et al. [53| ■ 
The results of Clementi and Hoffmann [5y| appear to be 
better reproduced by DFT in this case (Table |TTT|. A po- 
tential source of error in our case is the fact that we are 
using functionals based on the high-density limit of the 
HEG [3, 113, IIjI which may not be a good approximation 
for electrons farther from the nucleus. 



VI. CONCLUSION 

We developed an electronic structure method based 
on pair density functional theory that is simple to im- 
plement. We have shown that exact energy functionals 
can in principle be constructed in terms of the Hartrcc- 
Fock pair density, and that the representability question 
can be circumvented. Our formalism is closely related to 
the work of Levy on Hartree-Fock density matrix func- 
tional theory [35|. In applications approximations have 
to be developed to account for correlation. In our method 
correlation is accounted for by auxiliary pair-potentials 



between electrons in a Hartree-Fock setting, thus the 
method scales also as Hartree-Fock. The correlation 
approximations tested give quantitative agreement for 
atomic test cases with other methods. In DFT there are 
many possible options for correlation functionals which 
have been developed over the past few decades, often 
with specific chemical or physical situations in mind. The 
close agreement for the two simple PDFT correlation 
functionals developed and tested in this work indicate 
that there is potential in developing correlation function- 
als for PDFT as well. Since the pair density is a quantity 
depending on the coordinates of two particles it can be 
expected to give a better overall, more robust description 
of correlation than the ones used in one-body DFT. 
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APPENDIX 

In this appendix we show how to determine the relation 
between the constants for a correlation energy for the 
HEG of the form [1, d [H 



ec{p) = A + B\np. 



(25) 



The two-body functional form we assume for the corre- 
lation functional can be written as 



ec{nx) = A + 



B 

2N 



drr'hx{r, r')ln|r — r'| 



(26) 



Using the fact that 

n,(|r-r'|) = p2[5(|r-r'|)-l], 
we can rewrite Eq. (PS|) as 

e(nj.) = A + 2'KBp / r'^drh^{r)\n\r\. 

Scaling the coordinate r as 

X — kfr 



(27) 



(28) 



(29) 



results in 



kf = (3^V)^ 



2B 



ee(n,) ^ A+— x'dx[gix)-l]\n\x\ (30) 



B. 

"6 



IuStt p. 



Using the form oi g{x) [9[, the relation between the con- 
stants can then be shown to be 



A = A-— I x^dx[g{x) - l]ln|2;| (31) 



-Bln37r^ 



B = 6B. 
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